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Abstract. We investigate the effect of a density gradient on Chandrasekhar's dynamical friction formula based on the method 
of 2-body encounters in the local approximation. We apply these generalizations to the orbit evolution of satellite galaxies in 
Dark Matter haloes. We find from the analysis that the main influence occurs through a position-dependent maximum impact 
parameter in the Coulomb logarithm, which is determined by the local scale-length of the density distribution. We also show that 
for eccentric orbits the explicit dependence of the Coulomb logarithm on position yields significant differences for the standard 
homogeneous force. Including the velocity dependence of the Coulomb logarithm yields ambigous results. The orbital fits in 
the first few periods are further improved, but the deviations at later times are much larger. The additional force induced by the 
density gradient, the inhomogeneous force, is not antiparallel to the satellite motion and can exceed 10% of the homogeneous 
friction force in magnitude. However, due to the symmetry properties of the inhomogeneous force, there is a deformation and no 
secular effect on the orbit at the first order. Therefore the inhomogeneous force can be safely neglected for the orbital evolution 
of satellite galaxies. For the homogeneous force we compare numerical N-body calculations with semi-analytical orbits to 
determine quantitatively the accuracy of the generalized formulae of the Coulomb logarithm in the Chandresekhar approach. 
With the local scale-length as the maximum impact parameter we find a significant improvement of the orbital fits and a better 
interpretation of the quantitative value of the Coulomb logarithm. 
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1. Introduction 

Dynamical friction has a wide range of applications in stellar- 
dynamical systems. It appears as the first order term in the 
Fokker-Planck approximation (see e.g. Binney & Tremaine 
I1987I hereafter BT, Eq. (8-53)) competing with the second or- 
der diffusion terms to allow for a Maxwellian as the equilib- 
rium distribution. In the field of relaxation processes (dissolu- 
tion of star clusters, gravo-thermal collapse, mass segregation) 
the physical constraints and dominating aspects of dynamical 
friction are very different to those of the dynamical evolution 
of a single massive body in a sea of light background particles. 
In the latter case the friction force is independent of the mass 
of the background particles, and as long as the massive body 
is far from energy equipartition the diffusion can be neglected. 
Its main applications are the satellite galaxies in the dark mat- 
ter halo (DMH) of their parent galaxy, and super-massive black 
holes (SMBH) or compact star clusters in the central bulge re- 
gion of galaxies. Even in a collisional gas, dynamical friction 
occurs, but corrections due to pressure forces must be applied 
(Just et al. [T986l Just & Kegel [T990l Ostriker [T999l using mode 
analysis, and Sanchez-Salcedo & Brandenburg 1999 200Qus- 
ing numerical calculations). 



There are different approaches to analyse the parameter de- 
pendence and to calculate the magnitude of the friction force, 
most of them based on perturbative methods (2-body colli- 
sions: Chandrasekhar lT3551 Binnev [T977l Smtzer [19871 BT; 
mode analysis: Marochnik 1968 Kalnais 19721 Tremaine & 
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Weinberg 1984, theory of linear response: Colpi et al. 1999; 
statistical correlations in the fluctuating field: Maoz 1993 1. Due 
to its complicated structure, dynamical friction shows many 
different features. In a recent theoretical investigation Nelson 
& Tremaine (1999 1 connected many physical aspects of the 
fluctuating gravitational field on the basis of the linear re- 
sponse theory. For practical applications one basic question 
is whether dynamical friction acting on a massive object is 
dominated by local processes. There is a longstanding discus- 
sion on the contribution of global modes to dynamical friction, 
which are known to be excited and enhanced by the self-gravity 
of the perturbation (Weinberg [T986I [T989I fT9l)3i Hernquist & 
Weinberg 1989 ). These strong resonances may reduce the dissi- 
pation rate considerably (Kalnais [l972l Tremaine & Weinberg 
[T9841 Weinberg [T989l but see also Rauch & Tremaine ( fT996l 
for the opposite effect. For interacting galaxies with moder- 
ate mass ratios, the global deformation of the main galaxy and 
the strong tidal forces on the perturber can enhance the en- 
ergy and momentum loss by considerable factors, which can- 
not be described by the Chandrasekhar dynamical friction form 
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(Prugniel & Combes fl992l Leeuwin & Combes [T9971 and the 
discussion therein). The mass of the perturber should not ex- 
ceed ~ 5% of the mass of the host object in order to neglect 
both the global contribution to the dynamical friction force and 
the effect of deformation. 

This paper is concerned with the local contribution to dy- 
namical friction. We restrict our analysis to the low mass 
regime with M < 0.05 M^, where M and are the mass 
of the perturber and the host object, respectively. For clearness 
we formulate most equations in terms of the orbital motion of 
satellite galaxies in the DMH of the parent galaxy and give 
additionally the explicit application to a singular isothermal 
halo. This case will also be used for the numerical analysis. 
Nevertheless the results can be used for more general back- 
ground distributions with isotropic distribution functions and 
also for compact objects. As a second application we give the 
corresponding equations for massive Black Holes moving in 
galactic centres. 

The basic formula of Chandrasekhar relies on an infinite 
homogeneous background with an isotropic velocity distribu- 
tion function. For a Maxwellian velocity distribution function 
and a massive object with M 3> m it is given by 



^hom — 
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(BT Eqs. (7-14), (7-18)). Here M and % are mass and velocity 
of the massive body; m, po, and a are mass, mass density, and 
velocity dispersion of the ambient particles, respectively. In this 
formula, only 2-body interactions between the massive particle 
and the background are included, adopting a straight-line mo- 
tion in the unperturbed potential. Essentially all uncertainties of 
the relevant regime in phase space for these 2-body encounters 
are hidden in the Coulomb logarithm In A. Different perturba- 
tion methods lead to different interpretations of A: (1) 2-body 
encounters: A = & max /&mm with maximum and minimum 
impact parameter (BT) (2) Mode analysis: A = fc max //c m i n 
with maximum and minimum relevant wavenumber (Kalnajs 

A = 



I1972i (3) Theory of linear response: A = T max /r m i n with 
maximum and minimum time lag of the correlations (Nelson 
& Tremaine 1999 1. The most important correction to Eq.^is 
due to an anisotropic velocity distribution (Binnev H977l Statler 
1988] I199U . A numerical investigation of the anisotropy is 
discussed in another paper for satellite motion in flattened 
Dark Matter haloes (Penarrubia l20031 Penarrubia et al.|2004i. 
In different approaches the inhomogeneity of the background 
distribution is generally included (Kalnajs |l 9721 Tremaine & 
Weinberg Colpi et al. [T9^9l Nelson & Tremaine [T999 1 

but a useful scheme for applications is not yet available. An 
explicit formula for corrections due to the local density gradi- 
ent is given by Binney 



19771, who has shown that the addi- 
tional inhomogeneous force can be neglected for the dynamics 



of galaxies in galaxy clusters. Maoz ( 1993 1 and follow-up pa- 
pers ( Sanchez-Salcedo 1999 Dominguez-Teneiro & Gomez - 
Flechoso 1998 1 used the fluctuation-dissipation theory with ap- 



proximations similar to those presented here to derive the en- 
ergy loss due to dynamical friction. Del Popolo & Gambera 
(1999 1 and Del Popolo (2003 1 compute the corrections due to 
a density enhancement spherically symmetric with respect to 
the massive body, which may be used to estimate the shielding 
effect on the friction force. 

Numerical calculations suffered for a long time from res- 
olution problems (e.g. Lin & Tremaine 1983, White [1983 
Bontekoe & van Albada 1987 , Weinberg 1993 1 and, even in re- 
cent years, computations with particle numbers below 10 5 are 
frequently used (Cora et al. 1997] van den Bosch et al. lT999l 
Jiang & Binney 120001 Hashimoto et al. I2003i . Low resolu- 
tion raises two major problems for the determination of orbital 
decay times and quantitative numbers for the Coulomb loga- 
rithm In A. Due to noise, the numerically determined Coulomb 
factors A vary significantly for the same orbits. With differ- 
ent seeds and evolution times for relaxation of the unperturbed 
galaxy using 5,000 particles Cora et al. ([1997 1 found a varia- 
tion of fitted values for A by almost a factor of 3. Even in the 
20,000 particle runs the single orbits differ by up to a factor 
of 2 in the satellite energy. Secondly, increasing the resolution 
should lead to systematically larger values for the Coulomb 
logarithm, because the range of impact parameters, where the 
perturbation of the motion of the background particles is re- 
solved, increases. This is explicitely shown in Spinnato et al. 
(2003 1. Therefore even an average over many low resolution 
computations underestimates the friction force and leads to too 
large decay times for the satellite galaxies. With the new gener- 
ation of computers it is now possible to perform a larger num- 
ber of self-consistent calculations with particle numbers above 
10 6 and the corresponding high grid resolutions in particle- 
mesh codes (Velazquez & White 1999, Penarrubia et al. 2002 
Spinnato et al. 2003 1. In these calculations one can hope to sup- 
press noise sufficiently and to also resolve the small scale per- 
turbations to model the full range for the Coulomb logarithm. 
Bertin et al 



I2003t also used a large number of particles but 
only a low spatial resolution due to the restriction on low or- 
der spherical harmonics. For the excitation of resonant global 
modes in N-body calculations, one order of magnitude higher 
particle numbers may be necessary. But this is not our aim here. 

The standard usage of Chandrasekhar's formula works with 
an individual fitting of the constant Coulomb logarithm for 
each orbit. Fitting the Coulomb logarithm to a single orbit pro- 
vides in most cases an accurate description of a large part of 
the orbit. For this task, there is no need to improve the method. 
But the Coulomb logarithm varies systematically from orbit to 
orbit as discussed above and the orbital fits drift away for long- 



living satellites (Lin & Tremaine 1983 Jiang & Binney 2000 
Spinnato et al. 120031 Penarrubia et al. 2004 1. In order to get 
an improved friction formula, which can be used for a sample 
of orbits with different parameters, we apply local constraints 
yielding a systematic variation of the Coulomb logarithm with 
position (and also velocity). The distance to the galactic cen- 
tre as the maximum impact parameter for the orbital evolution 
of the Magellanic Clouds was already proposed by Tremaine 



( 1 976 1 and a position-dependent Coulomb logarithm was found 
numerically by Bontekoe & van Albada ( 1987 1 and Bertin et al. 
(2003 1. However the properties of a position-dependent In A 
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were never investigated systematically. We discuss also the ve- 
locity dependence of A by excluding slow encounters. 

A second effect of the local density gradient of the back- 
ground distribution is an additional 'inhomogeneous' force. 
Due to the different symmetry of the density gradient, the net 
effect of the perpendicular component of 2-body encounters 
does not vanish like that for the homogeneous part, but leads to 
an additional force not anti-parallel to the satellite motion (see 
also Binney 1977). We discuss the new Coulomb logarithm and 
the inhomogeneous friction term in detail. 

We calculate dynamical friction using the explicit deflec- 
tion of 2-body encounters, because this approach is not based 
on the diffusion limit in phase space. Therefore, in the case of 
a massive particle, the derivation also holds for large angle de- 
flections leading to a well-defined behaviour for small impact 
parameters and slow encounters. Neglecting the effect of the 
mean tidal field on the encounter yields the well known deflec- 
tion 



Al>M|| = 



2m Vo 



1 



A«M1 



M 
2Gm Vo 

2mbV§ 
GM 2 
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GM 

~V T 
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(4) 
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(BT, Eqs. (7-10a,b)) for a given impact parameter b and relative 
velocity Vo- Here a(Vo) corresponds to the impact parameter 
for a 90° deflection. From this acceleration we have to subtract 
the acceleration along the unperturbed orbit properly in order 
to avoid artificial results due to differences in the mean field 
approximation. 

We start in Sect.|2]with a detailed discussion of the collision 
rate of 2-body encounters and the structure of the Coulomb log- 
arithm with all parameter dependences. In Sect.[3]we compute 
the inhomogeneous dynamical friction for an isotropic distribu- 
tion function including a reanalysis of the homogeneous fric- 
tion force (corresponding to the standard friction) for a con- 
sistent comparison of both parts. In Sects. |4] and [5] we ap- 
ply the new formulae to typical satellite orbits in the Dark 
Matter halo of the Milky Way and discuss the effect of the im- 
proved Coulomb logarithm and of the inhomogeneous friction 
term. Sect. [6] contains a discussion of the relevance of using a 
position-dependent Coulomb logarithm and of the inhomoge- 
neous force for other applications (super massive black holes in 
galactic cores) and other aspects (circularization of the orbits). 

2. Encounter rates and the Coulomb logarithm 

In this section we discuss the encounter rate and perform the 
integration over impact parameters for the homogeneous and 
the inhomogeneous contribution to dynamical friction. General 



symmetries of the force and the properties of the Coulomb log- 
arithm are described in detail. 



2. 1. Encounter rate and mean field correction 

We use a local coordinate system (p, q, s) with s parallel to 
Vo and (p, q) the plane of impact parameters perpendicular to 
Vo (see Fig. I12> . The integration over s is substituted by the 
encounter rate n(r) \V~o\ with unperturbed background density 
n(r) and with r = rM + b using formally a 3-dimensional 
impact parameter in vectorial form, if necessary, given by b = 
(p, q, 0) = 6(cos 4>, sin <j>, 0). The rate of encounters for given 
relative velocity Vo and impact parameter b is then 



di/(6, Vb) = |Vo|n(r)/(«)dpdgd 3 V 



(6) 



the Dark Matter particles are counted when they are cross- 
ing the (p, g)-plane. The density gradient parallel to Vo cor- 
responds to an acceleration due to the mean field because the 
phase space density along (unperturbed) orbits is constant at 
equilibrium. The effect on the timing during the encounter 
can be accounted for only when tidal fields are taken into ac- 
count. Accelerations due to the mean field will not be included 
here. The density gradient perpendicular to Vb changes the en- 
counter rate directly. This effect of the local inhomogeneity is 
dominant and will be included. 

Since we are interested in the additional acceleration due 
to the grainy structure of the background distribution, we have 
to take the difference of the acceleration between the perturbed 
and the unperturbed orbit. Particles moving along the unper- 
turbed orbits with impact parameter b and relative velocity Vb 
correspond to a Tine' of constant density p(r) f (v)dsdpdqd 3 Vq 
along s, that is, parallel to Vq. By symmetry, only the parallel 
component to b contributes. We find 



i'mo = / Gp{r)f(v) 

2Gp(r)b 
b 2 



Ms 



> s/b 2 + s r 

f{v)dpdqd 3 V 



rdpdqd 3 V Q 



(7) 



for the acceleration along the unperturbed motion. This accel- 
eration must be subtracted from the acceleration due to 2-body 
encounters yielding an effective acceleration due to the grainy 
structure of the Dark Matter given by 



v di = J Av M dv(b, Vo) - vmo = #m|| + v M ± 



(8) 



Since the correction due to the unperturbed orbits is parallel to 
b, it affects only the perpendicular component. We find from 
Eqs.|3]andH 



wmii = 



«M1 



2Gp(r)f(v) Vq 
a(V ) Vo 

'2 Gp(r)f(v)b 
a(V Q ) 2 

2Gp(r)f(v)b 
b 2 



-i -l 



1 



a(Vo) 
b 2 



dpdqd'Vo , (9) 



a{Vof 
dpdqd 3 V 
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-2Gp(r)f(v)b 
b 2 



a(V ) 



dpdqd s V . 



For the integration over the impact parameters we expand the 
local density distribution up to the first order in the (p, q) -plane 



the perpendicular component) is not perpendicular to vm even 
for an isotropic distribution function. For vanishing vm and 
(10) isotropic velocity distribution function it is antiparallel to the 
density gradient. This is equivalent to an effective reduction of 
the enclosed mass M. (t*m) of the Dark Matter halo. 



p(r) w po + b ■ p± with r = r M 



(11) 2.3. Integration over the impact parameters 



Here po — p(r)\ TU and for the local density gradient we use 
the abbreviations 

Pr = Vp(r)|. rM = (p x ,p y ,p z ) = (p P ,p q ,Ps) and (12) 

P_L = Pr (Pr •e Vo )e Vo = (p p ,pg,0) (13) 

with coordinate systems (x, y, z) oriented to the velocity of the 
satellite with z axis parallel to %, and (p, q, s) with s axis 
parallel to Vb (both are locally centred at rM, see Fig. ll2t . The 
vector p± is the component of the density gradient perpendic- 
ular to Vq. In the Taylor expansion an additional variation of 
f(v) is neglected. So, for example, we do not account for some 
spatial variation of the velocity dispersion <r(r). 

2.2. Symmetries and general properties 

There are two important features of the velocity change due to 
2-body encounters: The symmetries and the mass dependence. 
Both are fundamentally different for the parallel and the per- 
pendicular component. 

Ai>m|| : In the limit of small angle deflections (b ^> a(Vb)) 
the parallel component is proportional to the mass M leading 
to the well known mass segregation. This also easily allows a 
separation from mean field effects. 

The parallel component is a function of b 2 only. Therefore, 
in a local Taylor expansion of the distribution function, only 
terms even in p and in q contribute to the integral over the 
impact parameters. These are the constant term, and diagonal 
terms of second derivatives, and so on. 

The parallel component is odd in Vq. Therefore the contri- 
bution vanishes for vanishing satellite velocity vm, if the distri- 
bution function f(v) is symmetric. 

At)M±: For small angles the perpendicular acceleration is 
proportional to M 2 leading to a stronger mass dependence. On 
the other hand, the acceleration falls off more steep with Vb 
leading to a dominating contribution at small Vb . 

The perpendicular component is odd in 6, therefore terms 
odd in p or q in the Taylor expansion of the distribution function 
contribute to the integral over the impact parameters, i.e. the 
linear term, 3 rd derivative terms, ... 

The perpendicular component is even in Vq. This has the 
effect of a net acceleration also for vanishing vm due to the 
asymmetric gravitational wake. 

Note that these symmetries are relative to the velocity Vb 
and not to the velocity % of the satellite through the Dark 
Matter. The directions of the accelerations relative to «m de- 
pend on the symmetries of the distribution function f(v). 
The zeroth order (parallel component), that is the standard 
Chandrasekhar friction, is only parallel to vm for an isotropic 
distribution function. The first order (the dominating term of 



Before integrating over the impact parameters we shall discuss 
the minimum and maximum impact parameter interpretation. 

Minimum impact parameter bo - . For compact objects the 
minimum impact parameter is irrelevant if bo <C a(Vb). In this 
case bo can be set to zero, but for a unified formulation we will 
use the Schwarzschild radius in the case of a Black hole. Since 
M 3> m, even in the case of large angle deflections, the veloc- 
ity change A«m is small and can be handled as a perturbation 
(i.e. under validity of using superposition of 2-body encoun- 
ters). For extended bodies, there is an additional effective bo 
determined by the typical radius of the satellite. The half-mass 
radius is a good choice with a correction factor of order 
unity, because the deflection angle is dominated by the force 
near the minimum distance and is determined by the enclosed 
mass. For a Plummer sphere the half-mass radius is given by 
rh = GM/ (4.6a 2 ) — 1.3r c with a central velocity disper- 
sion (To and core radius r c (the detailed calculation of White 
( 1976) corresponds to bo — 1.6r c ). An analytic calculation for 
extended bodies is given by Mulder ( 1983 1. For a King profile 
the half-mass radius is between 1/3 and 1/4.6 of GM/oq, again 
with central velocity dispersion oo of the satellite. Therefore 
we use for extended objects 



bo = r h 



GM 



(14) 



The latter approximation, which also corresponds to the half- 
mass radius of a singular isothermal sphere with outer cut-off, 
will be used for the investigation of the parameter dependence 
of the Coulomb logarithm. In case of a compact object, we use 
a characteristic radius for bo (the Schwarzschild radius r$ for a 
Black Hole) and define the corresponding virtual velocity dis- 
persion by ctq = GM /Abo leading for a BH to 



2GM GM 2 GM c 

bo = rs — — ^— = — — with a n 



Ao 2 



Abo 



8 



(15) 



in order to get a consistent description. 

Maximum impact parameter b\\ The allowed (and ef- 
fective) range of impact parameters is bounded by the geom- 
etry and by the relevant time-scales. The local approximation 
of the distribution function in density (and possibly in veloc- 
ity space also) breaks down for distances larger than the lo- 
cal scale-length L. Therefore the maximum impact parameter 
should be smaller than L, 



b < L = p/|Vp| = po/\Pr 



(16) 



By using this cut-off, we neglect more distant encounters and 
may underestimate dynamical friction. For a singular isother- 
mal sphere we have L = rjvi/2. We account for this uncertainty 
by using a fitting factor Q2 in the maximum impact parameter. 
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The straight-line approximation and the local approxima- 
tion require that the encounter time t co \\ = 2b/Vo is short 
compared to the local dynamical time-scale idym which is es- 
sentially determined by the local crossing time of the satellite 

idyn = L/ V M 

2b _ L 

T~ — ^coll ^ ^dyn — 

V v m 



(17) 



So the maximum impact parameter depends on V$. Using nor- 
malized velocities X = Vm/ (V^cr) and W = V~o/(V2cr), it is 
given by 



bi = Q 2 Lmm[l, —^-] 
2v M 



Q2L 



W 2 



AX 2 



W 2 



(18) 



In order to guarantee formally bi > fe also f° r verv small Vq 
we will use the analytic approximation 



b\ = bl + Q\U 



W 2 



AX 2 + W 2 



(19) 



Since again it is a priori not clear whether the time-scale re- 
striction does not exclude an important contribution of long- 
term encounters, we will also carry out computations using a 
maximum impact parameter given by 



2 T 2 



QfL 



(20) 



neglecting the time-scale argument. 



Fig. 1. The different approximations for the Coulomb loga- 
rithm as a function of normalized encounter velocity W = 
Vo/(^/2a). The full, dotted and dashed lines show ln(A2) 
(see Eq. 1281 for different normalized satellite velocities X = 
vm/{V2ct). The satellite parameters — 40, q s — 5 cor- 
respond to a mass of M — 5.6 ■ 10 9 M© at a distance of 
tm = 40 kpc in a Milky Way-like halo. For comparison the 
Coulomb logarithm ln(Ai) (see Eq.[3U| dashed-dotted line) is 
shown. The corresponding value for a Coulomb logarithm in- 
dependent of the encounter velocity ln(Ao) (see Eq. [33] thin 
line) is also plotted. We have chosen Qq = 0.8 and Q2 = 2.5 
from the numerical experiments and the corresponding value 
Qi = 1.0, which leads to the same friction force for X = 1. 



2.4. Homogeneous and inhomogeneous dynamical 
friction 

For clarity we split the dynamical friction v^f (Eq.[8} into the 
part induced by the homogeneous density term and that by the 
local density gradient (inhomogeneous friction) 



Vdf = Vhom + Vinh 



(21) 



By symmetry, one has that after integrating over the impact 
parameters, only the parallel component ■Um|| contributes to 
the homogeneous friction i>hom, and the perpendicular com- 
ponent vm± to the inhomogeneous friction ■Ui n h- For compar- 
ison we calculate the standard case (homogeneous part, paral- 
lel component) in the same way as the inhomogeneous term. 
Inserting into Eq. [5] the Taylor expansion of the density (Eq. 
II Q and integrating over impact parameters gives (the gradient 
term vanish by symmetry after integration over angle cf> from 
dpdq = d<j)bdb) 



with 



f(v) 



2G Po V 

a(V )V J h 



V 



-1 -1 



a(V ) 2 



irdb 2 d 3 V 



4ttG P o I ^a(V )HA)f(v)d 3 V 



(22) 



The inhomogeneous term arises from the leading contribution 
of the perpendicular acceleration and gives the dominant con- 
tribution to the acceleration perpendicular to Vb . The integra- 
tion over the angle <f> leads to 



/ (b ■ p±)bd(j) — irb 2 p± 
Jo 



(24) 



Using this result for the integration of Eq.^|over the impact 
parameters we find 



*inh = / nGp±f(v) 



1 



a(V ) 



d6 2 d 3 ^ 



2nG / p ± a(V y\n(A)f(v)d 3 V 



(25) 



where only the density gradient term of the Taylor expansion 
contributes. Note that p±, the density gradient perpendicular to 
Vb, is a function of Vb leading to a more complicated angular 
integral in velocity space. 

2.5. The Coulomb logarithm 

For the analysis of the parameter dependence of the Coulomb 
logarithm, we can reformulate Eq.[23]using two scaling param- 
eters, namely the size ratio q^ and the ratio of specific kinetic 
energies q s defined by 
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instead of the dependence on bo, b%, and the mass M. Since 
tTo is very large for compact objects (see Ea.ll5>. we get very 
small values for q s , which can in most cases be neglected. With 
Eq-El Eq.|5]becomes 



<V ) = ^ 



_2bo_ 

q s W 2 



(27) 



With Eqs.^Jand|23]we 8 et tne analytic approximation 



In A? 



In 



\ 



W± n 2 

A 'lS 



1 + Qbi 



w 2 



W 2 +iX 2 



1 + ^ 



Qhhi 



32X 2 

and alternatively with Ea.l20l 



for W ^ 



In Ai = In 



'1 + T9 S 2 (1 



1 + T^s 2 

for W ~> 



r\2 1 2 



W 



(28) 
(29) 

(30) 
(31) 



In order to isolate the effect of the dependence on the en- 
counter velocity W of the Coulomb logarithm, we also anal- 
yse In Ao independent of W by using the position-dependent 
maximum impact parameter b\ = QqL and an estimation 
ago of the typical impact parameter for a 90° deflection. With 
v 2 ^ inserted into Eq.[3]we find 

GM 2b 



(V 2 ) typ = 2a 2 



a-90 



2a 2 +v 2 M q s {l + X 2 



(32) 



which depends on tm and vm, but not on W. Inserting this into 
Eq.|23| we find 



lnA 



In- 



QoL 



In 



Qoqd 



00 



ln(Qo9d 
In 



Q L{2a 2 + v 2 M ) 
GM 



for extended objects 
for compact objects. 



(33) 

(34) 
(35) 



For the fits to the numerical orbits we also use a globally con- 
stant value (In A) . 

In the case of an isothermal halo we can substitute qa by 
using M.(ryi), the enclosed DMH mass, and the mass M of 
the satellite, where 



2 

Qm 



M(r M ) : 

M(r M ) 
M 



QsQd 



G 



For the isothermal halo the Coulomb logarithm becomes 



In A, 



In 



\ 



1 + ^k 2 



QU 



lV 2 



2 2 

2lmw 2 +4:X 2 



W± 2 
4 Hs 



(36) 
(37) 

(38) 



and similar for In Ai. In Fig.[^the velocity dependence of the 
Coulomb logarithm is shown for different satellite velocities X 
and typical parameters q^, q s (X = 1 corresponds to a circular 



orbit in an isothermal halo). We show the velocity-dependent 
Coulomb logarithm ln(A 2 ) for X = 0.1; 1.0; 3.0 (dotted, full, 
and dashed line). We have chosen Q2 = 2.5, which is the 
best fit value for typical orbits discussed below. For smaller 
X the cut-off at low relative velocities W is steeper, whereas 
for higher X it becomes flatter. The Coulomb logarithm ln( Ai ) 
(dashed-dotted line) with constant 6 max shows a shallower cut- 
off for W -> (oc W 4 for W 2 < 2/q s instead of oc W 6 
for W 2 <C 2/q s and W 2 <C 4A" 2 ), which occurs at a lower 
relative velocity. In general, the cut-off shifts the main contri- 
bution to the friction force to higher values of W. For the sim- 
pler approximation ln(Ai) we used Qi = 1.0 in order to get 
the same force for X = 1 after integration over W. Neglecting 
the velocity dependence by using ln(An) with the best fit value 
Qo = 0.9(thin line) leads to the corresponding constant value 
for comparison. The effect of the low velocity cut-off on the 
friction force is discussed in detail in Sect. [3] 

3. The dynamical friction force 

For the integration over velocity space we restrict the analy- 
sis to the isotropic velocity distribution functions f(v 2 ). An 
anisotropy plus a local density gradient would lead to very 
complicated angle integrations, which can hardly be performed 
analytically. All details of the integrations in spherical coordi- 
nates in d 3 Vb-space are given in the Appendix. Here we discuss 
the results of a Gaussian distribution 

f(v 2 ) = J_ exp(-^) (39) 

for the explicit integration, but the general results do not de- 
pend strongly on the special shape of f(v 2 ) and it can be done 
in the same way for other functions. 

The homogeneous friction vector Vhom is parallel to the 
satellite motion vm., whereas the inhomogeneous friction lies 
in the («m, Pr) -plane with angle ^ between i>m and p r . 
Therefore we decompose the inhomogeneous term into the par- 
allel and orthogonal components with respect to Vm. 



with 



(40) 



^inh — ^par i ^ort 

^par — ("^inh ' *^vm)^-vm 
■Wort = |*ort|e or t • 

It is convenient to use the normalized functions Ghom(^), 
G pai (X), and G OIt (X) 

2 



Ghom{X) 



cosh(2WA) 



1 poo 

= / ln(A) exp(-W 2 - AT 2 ) x 

JO 

smh(2WX)' 



dW 



G paT (X) — 



cosh(21¥X) 



2WX WX ' 

ln(A) exp(-W 2 - X 2 ) x 
dW 



V^Qdqs Jo 

sinh(2WX) 



G Q rt (X) 



1 



ainh(2WX) 



V^QdQs Jo 

AW 2 X 2 + 1 



2WX X 2 W 4 ' 

ln(A)cxp(-iy 2 - AT 2 ) x 
dW 



2WX 



cosh(2T^A:) 



X 2 W 4 



(41) 



(42) 



(43) 
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leading to 

*hom = -CchGhomPOevM ' ( 44 ) 

* P ar = -C C hG pat (X) cos(*)e VM , (45) 
*ort = -CchG or t(X)sin(*)e ort , (46) 

where the projection factors cos^) and sin(^E') of the density 
gradient (see Ea.l62> are separated. The scaling factor 1 / (qdq s ) 
(= l/q m for an isothermal halo) is included in the functions 
G par (A), Gort(X). The explicit derivation of these functions 
are given in the Appendix (see Eq.l77>. In a spherical halo the 
radial and tangential components v r , Vt of the inhomogeneous 
friction are 

v r = C ch [G ort (A)cos 2 (z) + G par (A)sin 2 (i)] e r , (47) 
«t = -Cch [G ort (A) - G par (A)] sin(i) cos(i)e t (48) 

with inclination i of vm with respect to et. 

3. 1. Homogeneous dynamical friction 

Here we shall discuss the effect of retaining the dependence of 
A on the encounter velocity W. Neglecting W in the logarith- 
mic term introduces an error of about 20% (see Fig. [3 in the 
friction force, which is a function of X. The velocity depen- 
dence on In A is twofold. The natural cut-off for large angle 
encounters due to the reduced efficiency is included in both ap- 
proximations, InAi and lnA2. In lnA2 additionally the slow 
encounter cut-off in the maximum impact parameter (see Eq. 
\l9i is included. The effect of these dependences can best be 
seen by comparing with the friction force using In Ao (Ea.l33> 
with velocity independent maximum impact parameter b\ and 
an average impact parameter ago for the 90° deflection. In Ao 
is still position and orbit dependent via L, a, and Vm- In Fig. [2] 
the force ratios using In Ai or In A2 with respect to In Ao are 
plotted as a function of satellite velocity X for different satel- 
lite parameters. The typical range of X for eccentric orbits is 
0.5 < X < 2. For an orientation the thin line gives the normal- 
ized friction force Ghom using In Ao. All other lines show the 
ratio of Ghom using InAi or lnA2 with respect to that using 
In Ao. We can learn three things from this plot. Firstly, the vari- 
ation is below 20% for all parameter sets. Therefore for a deter- 
mination of the magnitude of the dynamical friction force, the 
simple generalization In Ao in the Chandrasekhar formula will 
work well. Secondly, the systematic variation between apo- and 
peri-centre depends on the approximation and on the parame- 
ters. This means that for a detailed investigation of the evo- 
lution of satellite systems or of the circularization, it may be 
necessary to include the explicit velocity dependence in the 
Coulomb logarithm. Thirdly, there is a significant relative vari- 
ation of the strength of the forces along eccentric orbits. The 
ratio of the forces at apo- and peri-centre can differ by 30%, 
when using In Ao instead of In Ai or In A2. Since the amount 
of circularization of the orbit depends on the relative strength of 
the forces around apo- and peri-centre (see Sect. l5.1> . this may 
lead to essential differences in the orbital evolution of satellite 
galaxies. In Sects. |4]and[5]we gj ve some numerical examples 
to show that the new approximations of the Coulomb logarithm 
improve the orbital fits. 
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Fig. 2. Relative homogeneous friction forces as a func- 
tion of satellite velocity X for the different approxima- 
tions of the Coulomb logarithm. The thin line shows the 
friction force Ghom using In An for orientation. The other 
lines give the relative variation Ghom(Ai)/Ghom(An) and 
Ghom (A2) /Ghom (Ao), respectively, for two sets of parame- 
ters Od, 9s- The first case is for a satellite galaxy, whereas the 
second case corresponds to a point-like object (SMBH with 
M = 2 10 7 M Q ). 




1 1 — 1 1 

0.5 1 1.5 2 2.5 3 

X 



Fig. 3. The orthogonal component of the inhomogeneous 
force Gort(A) normalized to Gch 1S shown using the differ- 
ent Coulomb logarithms: InAi (Qi — 1.0, dashed lines) and 
In A 2 (Q2 = 2.5, full lines). For comparison the homogeneous 
force with velocity independent In A (Qo = 0.9, dotted lines) 
is also plotted. The functions are plotted for different parame- 
ters: gd = 40, q s = 5 thick lines; gd = 6, q s = 5 thin lines. 

3.2. Inhomogeneous dynamical friction 

The properties of the inhomogeneous force given in Eqs. H31 
and 146 1 are much more complicated. We first discuss G or t(X), 
the component orthogonal to the satellite velocity which pro- 
vides the inhomogeneous force in units of Gch, if the density 
gradient is perpendicular to the satellite motion (at peri- and 
apo-centre in a spherical halo). In Fig. [5] the contributions of 
G it using InAi and In A2 are compared to Ghom with In Ao 
(thick lines: g<j = 40, q s — 5, thin lines: g<j = 6, q s — 5). 
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Fig. 4. Ratio R le i(X) = G pai (X)/G OIt (X) of the parallel and 
orthogonal contribution to the inhomogeneous force (without 
the projection factors cos >J/ and sin for different parameters 
qd, <7s- 



We have chosen the same values Qo — 0.9, Qi = 1.0, and 
Q2 = 2.5 as before. With In Ai the dependence of G or t on X 
is relatively flat for low satellite velocities. In the most relevant 
regime 0.5 < X < 2, the force can exceed 15% of the ho- 
mogeneous force (thin lines in Fig. [3}- Using In A2 leads to a 
steep decline for increasing satellite velocity. For X w 0, the 
value is a factor of about 2 higher and for X w 1 it is a fac- 
tor of 2 smaller compared to the In Ai case. The magnitude is 
up to 10% of the homogeneous force for eccentricities not too 
large. The figure also shows that the parameter dependence is 
quite different to the homogeneous force. The dependence on 
the relative compactness parameter q s of In Ai and In A2 partly 
compensates for the reduction of due to a more extended 
satellite. 

The component parallel to the satellite motion described by 
G par is smaller than the orthogonal contribution (see Fig. 0}. 
The inhomogeneous force vector lies in the (»m, p r ) plane, 
and the angle r\ between v ln h and ■um is given by 



tan(?7) 



tan(>J/) 

Rrel(X) 



with 



Rrcl{X) 



Gp ar (A") 
Gort(A^) 



(49) 



For small X the force is approximately parallel to p r , because 
-Rrel(O) = 1, and for large X the force is approximately parallel 
to e or t with R ie i(X) oc X~ 2 . 

In a spherical halo the radial component points outwards 
along the whole orbit and the tangential component changes 
sign at apo- and peri-centre (see Eas.l47landl48>. Therefore due 
to the symmetry of the unperturbed orbits with respect to apo- 
and peri-centre there is no secular effect of the inhomogeneous 
force. The only effect is a slight deformation and enlargement 
of the orbital shape. 

Along orbits with high eccentricity, the inhomogeneous 
force can dominate over the homogeneous contribution around 
the apo-centre, where X <C 1. In this case there is a net accel- 
eration before reaching the apo-centre. The effect is an outward 
shift of both the apo- and peri-centre leading to a longer dura- 
tion of apo-centre passage relative to peri-centre passage com- 



pared to the orbit without inhomogeneous force. This may in- 
fluence the evolution of highly eccentric orbits significantly due 
to an enhancement of angular momentum loss. This scenario is 
very speculative, but illustrates the general properties of the in- 
homogeneous force. Other effects like the bending of the grav- 
itational wake along curved orbits due to the time lag may lead 
to a much stronger correction especially at apo-centre. This is 
not included in the instantaneous and quasistationary ansatz of 
the Chandrasekhar dynamical friction. 

4. Satellite motion in Dark Matter haloes 

In order to illustrate the effects that the different approxima- 
tions for the Coulomb logarithm induce on the satellite mo- 
tion in a DMH, we compare N-body calculations with semi- 
analytical computations of typical orbits. A comprehensive nu- 
merical study is under the way and will be presented in a 
follow-up paper. We choose parameters representative of those 
inferred for the Milky Way. We use a rigid satellite, represented 
by a point-mass, to avoid mass loss and deformation effects. 
This assumption may be somewhat academic, but for the anal- 
ysis of the influence of the different approximations for the 
Coulomb logarithm on the orbits it is a useful simplification. 
In this case the minimum impact parameter is determined by 
the numerical resolution and not the size of the satellite. 



4. 1. The N-body code 

We use Superb ox (Fellhauer et al. 2000) to evolve the galaxy- 
satellite system. SUPERB OX is a highly efficient particle-mesh- 
code with nested and comoving grids based on a leap-frog 
scheme, and has been already implemented in an extensive 



study of satellite disruption by Klessen & Kroupa ( 1998 1 and 
Penarrubia et al. d2002t . 

Our integration time step is 0.39 Myr which is about l/25th 
the dynamical time of the satellite. We have three resolution 
zones for the DMH, each with 64 3 grid-cells: (i) The inner grid 
covers the central 63 kpc in order to have a uniform cell-size of 
2. 1 kpc along the satellite orbits, (ii) The middle grid covers the 
whole galaxy with an extension 164 kpc giving a resolution of 
5.6 kpc per grid-cell, (iii) The outermost grid, which is common 
for all components, extends to 348 kpc and contains the local 
universe, at a resolution of 11.6 kpc. The satellite is modeled 
by an additional massive particle. 



4.2. The semi-analytical code 

We have developed a simple numerical algorithm to integrate 
the equations of motion 



Tm = -V$ g (r M ) + i>df ■ 



(50) 



where $ g is the galaxy potential and ■Udf the dynamical fric- 
tion acceleration acting on the satellite along the orbit. In this 
scheme, the DMH is represented by a fixed density profile 
Ph{j") (Ea.l51>. the local velocity dispersion is determined an- 
alytically from the Jeans equation. The satellite is represented 
by a point mass with mass M . With this kind of semi-analytical 
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run 


e 


r p 


r a 


(In A) 


Qo 


Q 2 






[kpc] 


[kpc] 








1 


0.0 


55.0 


55 


2.4 


0.6 


1.9 


2 


0.5 


18.0 


55 


2.2 


0.9 


2.5 


3 


0.8 


5.3 


55 


1.8 


0.9 


2.8 


mean 








2.0 


0.9 


2.5 



Table 1. The numerical experiments. The peri- and apo- 
galactica are r p and r a , respectively, and e = (r a — r p ) /(r a + 
r p ) is the initial orbital eccentricity. The best fit values for the 
different analytical approximations of the Coulomb logarithm 
are given in the last columns (taken from Fig. [8}. The mean 
values are over all eccentricities. 



ansatz, the offset of the DMH to the centre-of-mass (~ 1 kpc) 
is neglected. 



all with apo-centre at r a = 55 kpc and peri-centre ranging 
down to r p = 6.1 kpc. In order to check the effect of numeri- 
cal noise, we performed two additional runs for comparison. In 
detail we find: 1) Increasing the number of particles by a factor 
of 8 yields an increase of the best fit Coulomb logarithm of 8%. 
2) With a new realization and a different orientation of the or- 
bit with respect to the grids, deviations of the orbits are below 
1 kpc in the first 3.5 Gyr. Then the orbits run slightly out of 
phase, but show the same general behaviour. We conclude that 
the basic results discussed below are not significantly disturbed 
by noise. 

We start all runs at apo-centre ro = ryiit = 0) = 55 kpc 
with a tangential velocity v(t = 0) = eve with local circular 
velocity vc from the rotation curve at the satellite's initial dis- 
tance. We start with appropriate values for e to reach the desired 
initial eccentricities. 



4.3. The galaxy and satellite parameters 

In order to avoid the effects of the disc on the satellite evolu- 
tion, we employ a galaxy model formed purely by a DMH, with 
characteristics similar to those usually assumed for the Milky 
Way. The density profile, therefore, can be approximated by a 
quasi-isothermal sphere with total mass M^, core radius 7, and 
cut-off radius r cut given by 

pn = M h g exp(-r 2 /r 2 ut ) 
27r 3 / 2 r cu t r 2 + 7 2 ' 

where the normalization constant 

a = {1 - \/¥/?exp(/3 2 )[l - erfOS)]}- 1 = (52) 
1 + ^/3+(tt-2) / 9 2 +0(/3 3 ) 

having (3 — 7/V cu t. In our calculations we use Mh = 1.57 x 
10 12 M Q , r cut = 168 kpc, and 7 = 3.5 kpc. The initial distri- 
bution function in velocity space follows a Gaussian (Eq. I39l > 
with cr = 150 km/s and a cut-off at v = V e (r) = \J— 2$(r), 
the escape velocity of the halo, approximating dynamical equi- 
librium. The number of halo particles is N — 1.4 10 6 in order 
to get sufficient resolution. The large cut-off radius is necessary 
to reach approximately an isothermal profile in the radius range 
of 15 kpc< r <60 kpc. The local scale-length drops from I = 
r/2.0 at r — 10 kpc to I = r/2.2 at r = 55 kpc only, and the 
velocity dispersion from 145 km/s to 125 km/s. The satellite is 
modeled by a massive particle with mass M = 5.60 x 10 9 M Q . 
The effective minimum impact parameter when using a grid 
code is smaller then the grid cell size. For the semi-analytical 
analysis it seems reasonable to use for the minimum impact pa- 
rameter half the cell-size, that is 60 = 105 kpc (see Spinnato 
et al. l2003l and Just & Spurzem 2004 for a discussion). 

4.4. Numerical experiments 

We started an extensive numerical study to analyse the effect 
of the improved Coulomb logarithm on orbital evolution. Here 
we shall discuss our first results concerning the quality of the 
orbital fits and the different systematic behaviour with vary- 
ing orbital eccentricity. Table Q shows the orbital parameters, 



5. Results 

For a better understanding of the fitting results, we first anal- 
yse the effects that each analytical treatment of the Coulomb 
logarithm introduces in the orbital evolution of the satellite 
by only considering the dominant homogeneous friction force. 
We first discuss the general orbital behaviour of the orbits. 
Subsequently we analyse the differences that each dynamical 
friction approach induces on the orbits. In Sect. l5.3l we consider 
the effect of the inhomogeneous force on the orbit evolution. 

5.1. Orbital evolution 

The orbital decay due to dynamical friction can be best ob- 
served by plotting the galacto-centre distance evolution. For the 
intermediate eccentricity e = 0.5 (run 2) this is shown in Fig. 
[5] The numerical data are given as dotted lines. For comparison 
the semi-analytical orbits with a constant Coulomb logarithm 
(In A) = 2.2 (full line) and with In A 2 using Q 2 = 2.5 (dashed 
line), the respective best fit values, are plotted. 

For the first 4 revolutions (2.5 Gyr) all fits are reasonably 
accurate. The simulations and the fits differ systematically in 
the late phase, which is a well-known problem in numerical 
simulations (see e.g. Bertin et al. 2003 1. The enhanced orbital 
decay seems to be connected to the large mass of the satellite. 
A similar effect does not occur in numerical simulations with 
lower mass (e.g. Hashimoto et al. I2003i . One possible reason 
for this behaviour is the nonlinear interaction with the (live) 
central region of the DMH. 

Since our main goal here is to analyse the effect of the max- 
imum impact parameter, we restrict our investigation to the first 
4 orbits. In Fig.[6]these are shown in the orbital plane for run 2. 
There is a systematic delay in the precession rate of the N-body 
orbit compared to both semi-analytical approximations. The 
reason for this systematic difference is the motion and defor- 
mation of the live halo due to the offset of the centre-of-mass, 
which cannot be reproduced by an analytic approximation of 
the DMH. 

The Coulomb logarithm in the different approximations 
and therefore the dynamical friction force varies in a very 
different way along eccentric orbits. In Fig. [7] the different 
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Fig. 5. Evolution of the galacto-centre distance for run 2, the 
live satellite with orbit eccentricity e = 0.5. The numerical 
data (dotted line) are compared to the semi-analytical fits using 
a constant (In A) = 2.2 (full line), and lnA 2 with Q 2 = 2.5 
(dashed line). 
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Fig. 6. Same orbits as in Fig. [5] in the orbital plane, but only 
the first 2.5 Gyr until the fourth time apo-centre position are 
shown. 



Fig. 7. The Coulomb logarithms and the dynamical friction 
forces along the first part of the orbit of run 2 until reaching 
the peri-centre first time. Upper panel: Orbital variation of the 
Coulomb logarithms (In A) (dotted line), In Ao (dashed line), 
and In ki{W = X) (full line) using the best fit values for the 
run. Lower panel: The mean field force GM. (^m)/^m> tne no " 
mogeneous friction forces Whom for both Coulomb logarithms 
as in the upper panel, and the inhomogeneous force also 
with Q2 = 2.5 are plotted. All forces are normalized to the 
mean field force at apo-centre. 



value with W — X is plotted. The different shape arises from 
the additional velocity dependence discussed in Sect. 13.11 The 
lower panel gives the real forces along the orbit. For compari- 
son the mean field force GM {tm}/^ and the inhomogeneous 
force are also plotted. The most important effect of a position- 
dependent Coulomb logarithm is that the relative importance 
of the friction force around the apo- and peri-centre is different 
by up to a factor of 2 for the different approximations. 



5.2. Best fit values 

In order to quantify the accuracy of the analytical approaches 
and to fix the fitting parameter Q in the Coulomb logarithm, we 
use the mean square fit of the apo- and peri-centre distances by 
minimizing the parameter x defined by 



2k 

-Y 

2k f-' 



(n 



(53) 



forces along the orbit of run 2 are shown until reaching the 
first peri-centre. The upper panel shows the Coulomb loga- 
rithms (In A) and In An. The decreasing value of In An with 
decreasing galacto-centric distance represents the position de- 
pendence of the local scale-length L. In the case of In A2 the 



7"j being the galacto-centric distances of the peri- and apo- 
galactica in the semi-analytical calculations. The subindex n 
denotes the corresponding values of the N-body code. The sum 
is over a given number of orbits k. The temporal off-sets Ati 
are also taken into account weighted with the velocity disper- 
sion at apo-centre. 
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Fig. 8. Fitting factor x as a function of the free parameters 
(In A), Q , and Q 2 , respectively. The different rows show the 
result for different eccentricities. The measure of \ is done for 
two numbers of orbits, k — 3, 4 (full and dotted lines, respec- 
tively). The last row gives the average over all eccentricities. 
We find well-defined minima to determine the best fit values 
listed in Tabled 
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Fig. 9. Secular evolution of specific energy, angular momen- 
tum, and of the eccentricity e for the orbit of run 2. The full 
curve shows the case of constant Coulomb logarithm (In A), 
the dashed line that with In A2, and the dotted line the N-body 
simulation. The values of e are given at peri-centre (triangles 
for the constant Coulomb logarithm and squares for the vary- 
ing Coulomb logarithm). The dotted line connects the numeri- 
cal values. 



The resulting x-curves as a function of adopted Q-values 
are plotted in Fig. [8] The upper rows give the values for the 
different eccentricities and the last row is the average over all 
available eccentricities given in Tab.[Q If the analytical formu- 
lae are perfectly correct, the best fit Q-values should be equal to 
unity and should be independent of eccentricity. In general the 
\- values show well-defined minima showing that the Coulomb 
logarithm is a well defined quantity. There is no big difference 
when using 3 or 4 orbits for the fit. We find an increasing qual- 
ity of the fit when using (In A), Qq, and Q 2 - In order to under- 
stand the systematic variation of the Q-values, we look first at 
the constant values (In A) = 2.4, 2.2, 1.8 for the different ec- 
centricities. The best fit value of (In A) decreases with increas- 
ing eccentricity as expected, since the average galacto-centric 
distance decreases. Looking at Qo, we find that the systematic 
variation with eccentricity is much smaller. From Qo ~ 1, we 
conclude that the maximum impact parameter is very similar to 
the local scale-length L. For Q 2 the trend changes and it seems 
that the radial variation is overestimated. The large value of 
Q 2 ~ 2.5 shows that an essential part of the friction force is cut 
off by using A2 with the time-scale argument. Slow encounters 
contribute essentially to the friction force. On the other hand, 
we find with I11A2 an even better fit than with InAo, which 
is already better than using the constant Coulomb logarithm 
(In A). Due to the complicated structure of In A2, conclusive 
results for the correct fitting formula (yielding Q 2 ~ 1) can be 
found with the full numerical investigation only. 



One important goal is to improve the parameter dependence 
of the Coulomb logarithm in order to use it for larger sets of 
orbits for a statistical analysis of the distribution of satellite 
galaxies. The variation of the best fit values of Q over the ec- 
centricity range yields a smaller variation of the corresponding 
Coulomb logarithms than in the case of using (In A). 

The evolution of specific energy E and angular momentum 
L due to drag forces depends on the eccentricity of the orbits. 
Using E = fM^df and L = ryiv^ cosi = Lv^vm with orbit 
inclination i we can write the ratio of energy and angular mo- 
mentum loss due to any drag force parallel to i>m in the form 

i = (54) 

which holds for all points along the orbit independent of the 
shape of the orbit. The energy loss is strongly enhanced at the 
peri-centre, where the satellite velocity is greater. The result 
for run 2 with e = 0.5 is shown in Fig. [9] The binding en- 
ergy increases in steps at the peri-centres, whereas the angular 
momentum decreases quite smoothly. The fit of angular mo- 
mentum is improved when using lnA 2 instead of (In A). The 
evolution of the shape of the orbit, measured by the eccentricity 
e, depends on the averages of energy and angular momentum 
change over the whole revolution from apo- to apo-centre. If 
the energy loss dominates, the orbit will be circularized, and 
if angular momentum loss dominates, the eccentricity will in- 
crease. The orbit with a constant Coulomb logarithm shows a 
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Fig. 10. Ratio of inhomogeneous to homogeneous force 
^inh/^hom (analytical from Eas. 1441-1461 along the orbit (full 
line) of run 2. For a better understanding we also plot radius 
Tm/tq (dotted line) and normalized satellite velocity X(t) of 
run 2 (dashed line). 

significant circularization, which is not observed in the N-body 
results (triangles compared to dotted line in the lower panel 
of Fig. |9}. The varying Coulomb logarithm essentially solves 
this problem leading to constant shape of the orbit (squares in 
the figure). The enhanced energy and angular momentum loss 
of the N-body calculation at later times will not be discussed 
here. 

5.3. Inhomogeneous dynamical friction 

Despite the problem of determining the maximum impact pa- 
rameter, the inhomogeneity leads to the additional acceleration 
Vmh, which has, due to the different parameter dependence and 
direction relative to satellite motion (see Sect. 0, a different 
physical effect on the orbits. 

Since the different contributions to the force on the satellite 
cannot be seperated in the numerical calculations, we discuss 
first the semi-analytical forces expected along a typical orbit. 
InFig.[lO| we plot the ratio of the inhomogeneous and homo- 
geneous force Winh/^hom along the orbit (run 2 with eccen- 
tricity e = 0.5). The maxima occur at the apo-galactica. The 
small bumps around the peri-galactica demonstrate the angle 
dependence of the inhomogeneous force being maximal, if the 
density gradient is perpendicular to the satellite motion. The 
magnitude of the forces are shown in Fig.0 The mean value of 
the ratio slightly increases as the satellite sinks into the inner 
region of the galaxy, but it remains below 20%. If we compare 
the semi-analytical orbits with and without the inhomogeneous 
force (see Fig. II It . the orbit correction is smaller than 1 kpc 
and dominated by a deformation of the orbit. The main effects 



r„ (horn) r„(inh) 




-0.5 - 

II. 

2 4 6 

t (Gyr) 

Fig. 11. Effects of including the inhomogeneous friction term 
Vj n h on the orbital evolution. Since the effects are small we 
show the galacto-centric distance difference of both orbits ArM 
(full line) and the difference in in 3-D space \rm{hom) — 
rm{inh)\ (dashed line) both in kpc. The dotted line shows the 
normalized position rM/Vo of the satellite. 



are a slight shift of the orbit to larger radii and a slow secular 
prolongation of the orbital time. For an understanding of this 
small effect we can look at the decomposition of the inhomo- 
geneous force into the parallel and orthogonal component (Eqs. 
I45land l46l i. The parallel component can be considered as a cor- 
rection to the homogeneous friction force, which changes sign 
when crossing the apo- and peri-galacticon. Therefore the main 
effect is a slight outside shift of the orbit. This can also be the 
reason for the slightly increasing orbital time of the satellite. 
There is no secular energy and momentum loss to first order. 
The orthogonal force leads to an additional bending of the or- 
bits. It points outwards and reduces the curvature of the orbit. 
The decomposition into the radial and tangential component 
(Eas.l47landl48> shows that the inhomogeneous force reduces 
the effect of the mean field force. The radial part works like 
a correction of the mean field by a few percent leading to a 
slightly different orbital shape with no secular part. 

Orbit deviations of this order of magnitude cannot be ex- 
tracted from the numerical calculations, because even with 
more than 10 6 particles the force fluctuates due to the noise 
in the halo particle distribution. From the small effect on the 
orbits we conclude that the inhomogeneous force can safely be 
neglected for the analysis of the kinematics of satellite galax- 
ies even if the magnitude of the inhomogeneous force is con- 
siderable. This confirms the conclusion of Binney dl977 i) for 
satellite galaxies. 
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6. Discussion 

We have recomputed the dynamical friction force on small 
satellite galaxies moving in the DMH of their parent galax- 
ies. We have extended the local approach of Chandrasekhar by 
adding local constraints for the determination of the friction 
force. We analysed the effect of the local density gradient of 
the DMH on the standard (homogeneous) friction force and the 
properties of the additional inhomogeneous force term. 

6. 1. The Coulomb logarithm 

We have shown that the local scale length L is an appropriate 
maximum impact parameter in the Coulomb logarithm In An. 
A numerical test for an isothermal DMH shows a significant 
improvement of orbits with different eccentricities compared to 
the standard case of a constant (In A) . This fits to the finding of 
Hashimoto at al. ((2003 ), who also used a position-dependent 
Coulomb logarithm. The numerical results of high-resolution 
simulations of Spinnato et al. ( 2003 1 are also in agreement with 
L as the maximum impact parameter. For a detailed analysis 
see Just & Spurzem (2004i. The improved Coulomb logarithm 
is also in agreement with the detailed investigation of Bontekoe 
& van Albada ( 1987 > for a n=3 poly trope with radius i? Q ut- 
They found for the range of central distances, where the density 
profile is nearly exponential (corresponding to a constant scale 
length L w 0.1i? ou t), constant values for In A, which are of the 
order of unity for a numerical resolution of ~ 0.04i? out . 

In general the numerical determination of In A strongly 
depends on the numerical method and resolution. In Cora et 
al. (1997) and other work a high level of noise and a wide 
spread in the numerical determination of the Coulomb loga- 
rithm are present. An extreme case is presented in Jiang & 
Binney (2000), who found In A = 8.5 from orbital fits to the 
Sagittarius dwarf, which corresponds to the unphysical ratio 
of fomax/frmin = 5,000. They used high particle numbers but 
an expansion into low order spherical harmonics with respect 
to the halo centre, which strongly suppresses the local gravita- 
tional wake, and artificially enhances the excitation of global 
modes by 'ghost satellites' due to the symmetry. Cartesian 
meshes are best adapted to a good representation of a local per- 
turbation with well-defined resolution - here global modes may 
be underestimated by the lack of general symmetry. 

The generalization to a velocity-dependent Coulomb loga- 
rithm by an additional cut-off for near but slow encounters us- 
ing In A2 gives an ambiguous result. On the one hand we find 
an additional improvement of the orbit fits for the first reolu- 
tions. But on the other hand in the late phase the orbital decay 
in the numerical calculations is even stronger than with a con- 
stant Coulomb logarithm. The reason may be a deformation 
of the central halo region due to the relative large mass ratio 
q m . The linear decay when using In A2 is similar to the find- 
ings of Hashimoto et al. (2003). Additionally, the fitting factor 
Q2 ~ 2.5 shows that by cutting off slow encounters a signifi- 
cant contribution to the friction force is lost. The rescaling by 
using Q2 > 1 may change the parameter dependence essen- 
tially. Here further investigations are necessary for conclusive 
results. 



6.2. Circularization 

One of the most important effects of a position-dependent In An 
is the reduction of the circularization of the orbit. The variation 
of In A along the orbit is a sensitive quantity for the evolution 
of the orbital shape, since the effect is significant even in the 
case where the dynamical friction force is changed by about 
15% along the orbit. The circularization depends on the sys- 
tematic parameter dependence of the friction force along the or- 
bit, because the energy loss relative to angular momentum loss 
is strongly enhanced around peri-centre due to the higher satel- 
lite velocity (see Ea.l54>. The general effect of a reduced fric- 
tion force around peri-centre is a reduction of the energy loss 
leading to a lower circularization. The effect of the position- 
dependent Coulomb logarithm on the circularization of the or- 
bits is in any case important for eccentric orbits. The general 
finding of a constant or increasing eccentricity in N-body sim- 
ulations can be understood in this way. 

6.3. Inhomogeneous dynamical friction 

The inhomogeneous term describes the effect of the asymme- 
try of the gravitational wake with respect to the orbit. For the 
calculation of the inhomogeneous force it is essential to use 
the velocity-dependent lnAi or lnA2 in order to have a ro- 
bust and consistent cut-off at low impact velocities, where the 
dominant contribution arises (see Appendix for a comparison 
with Binney ( 1977), who used a lower cut-off velocity). The 
inhomogeneous force is for typical parameters of the order 
of 10% relative to the homogeneous force. In highly eccen- 
tric orbits around apo-centre and for higher satellite masses 
(~ 10 10 Af Q ) it becomes significantly stronger. The force ratio 
of inhomogeneous and homogeneous terms scales essentially 
with ~ M/M.(ryi) (the ratio of satellite mass and enclosed 
mass of the DMH). This is very similar to the ratio of the or- 
bital to the friction time-scale. 

The properties of the inhomogeneous term are very differ- 
ent to the standard friction force. For an isotropic velocity dis- 
tribution it is an even function of the satellite velocity leading to 
a non-vanishing force at vm = and to a non-secular behaviour 
around apo- and peri-centre. The inhomogeneous force is less 
inclined to the satellite motion than the density gradient. For 
vanishing satellite motion it points anti-parallel to the density 
gradient. The parallel component gives an acceleration during 
outward motion and a deceleration during inward motion. The 
net force points outwards of the orbit leading to a less eccentric 
orbit. 

In our test orbits, the effect of the inhomogeneous force is 
very small leading to an orbital offset of the order of 1 kpc over 
a few Gyr. Since the magnitude of the inhomogeneous force is 
so small and their is no secular effect to first order, the correc- 
tions to the friction force must be compared to other shortcom- 
ings of the dominating homogeneous force. In highly eccen- 
tric orbits the time lag of the perturbation leads to a stretching 
or squashing at peri-centre and a bending at apo-centre of the 
gravitational wake due to the accelerated unperturbed motion 
of the satellite. The corresponding corrections may affect the 
orbital evolution much more than the inhomogeneous force. 
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6.4. SMBHs in galactic centres 

For compact objects an additional velocity dependence of the 
Coulomb logarithm via ago, the typical impact parameter for 
a 90° deflection, appears in InAo. For an application to the 
orbital decay of super-massive Black holes (SMBH) in cuspy 
galaxy centres, a significant correction to the decay time-scales 
is required. A quantitative analysis of the motion of SMBHs in 



galactic centres is given in Just & Spurzem (2004i. 
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Appendix: Integration over velocity space 

Coordinate systems 

We discuss the simplest case of an isotropic velocity distribu- 
tion function f(v 2 ). Since the angle integration in v-space is 
very complicated, for the inhomogeneous term at least, we use 
spherical coordinates in Vq -space oriented at vm- This leads to 



d 3 Vb 

„,2 



= -V 2 d<pdpdV , 



cost 



= V Q (y/l - p 2 cos ip, \J\ - p 2 sin<^,/i) , 



P± = 



\V - 

Py 



dm 



Vq 2 + 2pV v M + v 



2 

M ' 



(55) 
(56) 
(57) 



cos 



tp \J 1 - p? [p x cos ipyl - p 2 + 



p y sin <pyjl - p 2 + p z p] , 
- sin^vl — P 2 [px cos ipyl—pfi + 
pysmp^jl - p 2 + p z p], 



Pz - p[Px COS ip-\/ 1 - fl 2 + 

p y sin tpy/l - p 2 + p z fj]\ 
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(58) 



(V2^r) 3 
with u = 2WX 



exp(-X A - W z ) exp(-up) 
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(67) 



In the standard case the angle integrals are greatly simplified for the explicit integration, but the general results do not de- 

by using the general relation pen d strongly on the special shape of f(v 2 ). For the homoge- 

H H neous term (Ea.l63> we need the integral 

h(V )ev =—H(V ) = H(\v-v M \) (59) 

dVo dv M r L lcosh(u) smh(u) 



with H'(V Q ) = h(V ) 



and then integrating over w-space. But this method is not pos- 
sible for the inhomogeneous term, because the dependence of 



27T/Ltexp( — M/Lt)d/i = — 4-7T 



-47TM 



(68) 



To f01 ' " 1 



n ± on V is more complicated. Therefore we integrate also for and for the i nhomoge neous components (Eqs. El and |65} the 



the homogeneous term directly over the angles. 

Integration overip 

Since v 2 is independent of ip, we can integrate Eqs. 1221 and!25l 
over ip without specifying the distribution function f(v 2 ). We 
use 

e Vo d</3 = 27r(0, 0, ft) = 2npe VM , (60) 



integrals 



J 2w[l — p 2 ] exp(— up)dp — 

( cosh(it) sinh(u) 
\ u 2 u 3 



8tt 



(69) 



8tt 

T 



To J 



for u < 1 
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p ± dip = 7T ([1 + ft 2 ]p x , [1 + ft 2 ]p y , 2[1 - p 2 ]p z 
vj- ([1 + p 2 } sin(*)e ort + 2[1 - p 2 } cos(*)e Vlk 



7r[l + /i 2 ] cxp(— up)dp 



1(61) 



-47T 



cosh(u) 



1 + u 2 



sinh(w) 



(70) 



with L = pq/IPj-I an d the decomposition of the density gradi- 
ent with respect to vm 



— 1 + — for u < 1 



Pr = |Pr| [cos(<J>)e VM +sin(*)e or t] , 



where for completeness the Taylor expansions for small veloc- 
(62) ities are also given. 

Inserting the results of the angle integration into Eq.[64] we 



where e or t is the unit vector in the direction of the orthogo- 
nal component of p r and <J/ the angle between v M and p r . find for the homogeneous term 
We split the inhomogeneous term in the same way into the par- 4irG 2 Mpo % 

allel and orthogonal component (see Eq. I40t . We find for the 
homogeneous acceleration term and the components of the in- 
homogeneous acceleration 



^hom 



7TO"' 
2 v2\ 



In A x 



^hom " 



AnGp Q e VM x 

2na{V ) HA)f(v 2 )pdpV 2 dV , 



exp(-W z - X 



-Cche v 



(63) 



vm Jo 

cosh(2WX) sinh(2WX)\ 
WX 2W 2 X 2 ) 

\n{X)g hom (X, W)dW. 



dW 



(71) 



Here we used Cch from Eq.|2]and the function 



-2irGp Q 

z 



COs(^)e* 



g hom (X,W) = -^exp(-W 2 -X 2 ) 



cosh(2WX) 



WX 



2ira{V ) 2 ln(A)/(^ 2 )(l - p 2 )dpV 2 dV Q , (64) 



siah(2WX) \ 



(72) 



-2irGp . , . 

■ sm(^)e ort x 



L 



ira(V ) 2 HA)f(v 2 )(l + p 2 )dpV 2 dV . 



2W 2 X 2 J 

(In A is a function of V 2 , too, see Ea.l23>. 

For the inhomogeneous terms we get from Eq.[64]the par 
(65) allel component 

-G 3 M 2 ^p r 



The integral over ft 



<j±X 2 



• cos( i S)e^ 



Since v 2 of the Dark Matter particles is a function of ft, we 
must specify the distribution function before going on. We use 
a Gaussian 

o,2 



ln(A) exp(-W 2 - AT 2 ) x 

sinh(2WX)\ dW 



1 



(V2¥o-) 3 



exp(~^) 



(66) 



cosh(2WX) 
= -C C hCos(*)e V] 



2WX J W 4 
ln(A)g pai (X,W)dW , (73) 
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Fig. 13. Here we compare the contribution to the homogeneous 
dynamical friction as a function of encounter velocity W for 
different satellite velocities X = 0.7,1.5 using the different 
approximations for the Coulomb logarithm. We chose values 
Qo = 0.9 (corresponding to the standard dynamical friction; 
thin lines), Qi = 1.0 (normal lines), and Q2 — 2.5 (thick lines) 
leading all to the similar homogeneous friction for X = 1. The 
system parameters are q d = 40 and q s = 5. 

and from Eq.[65]the orthogonal component 

Vort = ^4^2 sm(*)e ort x 

ln(A) exp(-W 2 - X 2 ) x 

AW 2 X 2 +1 \ dW 

2WX ™H2WX)-co sH 2WX)) w 

poo 

= -Cch sin(tf )e ort / ln(A) 5ort (X, W)dW , (74) 
Jo 

where we used analogous functions 



9ort(X,W) 

h(X, W) 



1 g hom (X,W) 
q d q s W 3 X ' 

1 2h(X,W)-g hom (X,W) 
q d q s 2W 3 X 

4= exp(-W 2 - X 2 ) sivh{2WX) . 

\/7T 



(75) 
with (76) 



Now we can discuss the effect of the different approxima- 
tions for the Coulomb logarithm (Eqs. [28] |30| I33i . For the 
homogeneous force (arising from ln(A)g} iom (X, W)) the ef- 
fect of the different Coulomb logarithms on the contributions 
as a function of W is shown in Fig. [O] The main effect of 
a varying Coulomb logarithm (with Qi or Q2) is a shift of 
the dominant contribution to higher encounter velocities. But 
the overall deformation remains small and the dependence on 
the satellite motion is weak. Therefore one can use for the 
homogeneous friction force the approximation of a velocity- 
independent Coulomb logarithm In Ao. The relative variation 
with X and dependence on the system parameters qd,q s is 
shown in Fig. EJ 

For the inhomogeneous force the effect of the Coulomb 
logarithm on both components (g pal (X, W) and g OI t(X, W)) 



X=1 ■ 
X=0.5 ■ 
X.2 ■ 



Fig. 14. Here the function <7d9sffort(A r , W) (see Ea.l76> for the 
orthogonal component of the inhomogeneous dynamical fric- 
tion as a function of encounter velocity W for different satellite 
velocities X = 0.5, 1.0, 2.0 is shown. It must be multiplied by 
the Coulomb logarithm, before integrating over W. The effect 
of using different approximations for the Coulomb logarithm 
is essential, because the dominating contribution comes from 
small W, where g OIt (X, W) has a singularity. For imagination 
the Coulomb logarithms with qa — 40, q s — 5, and X = 1 
from Fig. [2 are plotted again (with Qo = 0.9, Q\ = 1.0, 
Q2 = 2.5 thin, normal, thick dotted line, respectively). 



X=0.7 A, 

A 2 

X.1.5 A. 




Fig. 15. The function In(A)<7 or t(X, W) for the orthogonal 
component of the inhomogeneous dynamical friction as a func- 
tion of encounter velocity W for different satellite velocities 
X = 0.7, 1.5 is shown. The difference from using Ai (thin 
lines) or A2 (thick lines) is large and depends also on X and on 
the system parameters qa, q s . The same values as in Fig.ll3lare 
used. 

is similar. In Fig.[3]the parameter dependence of the Coulomb 
logarithm on the cut-off is demonstrated. The results on the 
contributions to the orthogonal component are shown in Fig. 
[2]using the same parameters as in Fig.[0] The main effect of 
a varying Coulomb logarithm (with Qi or Q2) is a well defined 
lower cut-off consistent with the determination of g OTt (X, W) 
itself. The different approximations give significantly differing 
results. Therefore we will discuss both cases in detail. Using 
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Fig. 16. Here the ratio g pax (X, W)/g OT t(X, W) of the contri- 
butions to the parallel and orthogonal component of the inho- 
mogeneous dynamical friction as a function of encounter ve- 
locity W for different satellite velocities X = 0.3, 1.0, 2.0, 3.0 
is shown. 

the approximation of a velocity independent Coulomb loga- 
rithm In Ao can only be used together with a somewhat arbi- 
trary cuttoff. This is not appropriate for applications due to the 
strong sensitivity of the resultant force on the choice of the cut- 
off value (see below). 

In Fig. ^] we show the contribution to the parallel force as 
a function of W relative to that for the orthogonal component. 
The parallel force is smaller for higher satellite velocities X, 
because the contribution from higher encounter velocities W 
is reduced with increasing X. Since the Coulomb logarithm 
determines mainly the contribution from small W, the effect 
on the orthogonal and the parallel components are similar . 

Integration over V 

We introduce now also the integrals of the normalized functions 
defined in Eqs.|72]|Z5]|76]by 



om,par,ort 



(X) = / ln(A) ffh , 



om,par,ort 



{X, W)dW (77) 



to get a convenient form for a discussion of the structure of the 
different contributions to dynamical friction (see Eqs. |^] P31 
I46l l. With a velocity-dependent Coulomb logarithm, the inte- 
gration over W must be done numerically. Here we give the 
derivation of the standard friction formula and a short discus- 
sion of the behaviour of the inhomogeneous terms at low values 
of W. The explicit form of the integrals for the different con- 
tributions are given in Sect. [3] (Eqs. |^ and 1431. Further 
discussion of their properties is also given there. 

Chandrasekhar's friction formula 

With the approximation of a Coulomb logarithm In An inde- 
pendent of W (Eq. 1331 we find the standard Chandrasekhar 
friction formula by integrating gh am (X, W) (Eq.l72l by parts 



g hom (X, W)dW 



exp( 



-W 2 -X 2 ) x 



co^WX) sinh(2WX)\ 
WX 2W 2 X 2 J 



dW 



V^x 
i 

+ ~x 



exp(-W 2 - X 2 ) 



2 sinh(2WX) 



2WX 



i f°° 

- I cxp(-T^ 2 - X 2 )smh{2WX)dW 
x Jo 

cx P (-x 2 ) + ^ J x cx P (-y 2 )dr 



l 

x 2 



erf(X)-^exp(-X 2 ) 

\/7T 



Inserting this into Eq.[^we nn d the result (Eq.^ 



Vhc 



-C C h In An 



X 2 



erf(X)-— C xp(-X 2 ) 



(78) 



(79) 



Since the maximum impact parameter depends on the local 
scale-length, the Coulomb logarithm depends on the position 
of the satellite. Neglecting this dependence also by using the 
global value (In A) leads to the standard formula. 

Inhomogenous force and the lower cut-off 

The functions g pllv (X, W) and g ort (X, W) show a singularity 
ex W~ 2 for vanishing W (Eqs. I76> . Since the velocity- 
dependent Coulomb logarithms with Qi and Q2 vansih for 
small W, i.e. lnAi cx W 4 and lnA2 oc W 6 , there is a nat- 
ural effective lower cut-off for the integration over W. But the 
main contribution arises from low W just above this cut-off. 
This leads to differences of a factor of two at X = 0, when 
using the corresponding best fit values Qi = 0.9 and Q2 =2.5 
for the orbit. 

If we want to use a velocity-independent Coulomb loga- 
rithm In An, we must introduce a lower cut-off Wo- With a 



similar approximation as used by Binney (1977|Eq. (7)), we 
find for X <C 1 the corresponding formula 



ah(X) ~ Cch 



8 1 (In A) 
3a/7T <7d<7s Wq 



exp(-X 2 ) . 



(80) 



Inserting the parameters at pericentre of our standard orbit = 
6, q s = 5, (In A) = 2.1, and G ort (A" = 0) = 0.5, we find the 
cut-off value Wo — 0.21. We may compare this value with the 
relative velocity Wu for a 90-degree deflection at the half-mass 
radius r^, which is 



W h = 



GM 
2<J 2 r h 



0.63 



(81) 



for the parameters above. Using Wh would decrease the inho- 
mogeneous friction force by a factor of three. 



